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Abstract 

We propose a multiscale method for elliptic problems on complex domains, e.g. domains 
with cracks or complicated boundary. For local singularities this paper also offers a discrete 
alternative to enrichment techniques such as XFEM. We construct corrected coarse test and 
trail spaces which takes the fine scale features of the computational domain into account. The 
corrections only need to be computed in regions surrounding fine scale geometric features. 
We achieve linear convergence rate in energy norm for the multiscale solution. Moreover, 
the conditioning of the resulting matrices is not affected by the way the domain boundary 
cuts the coarse elements in the background mesh. The analytical findings are verified in a 
series of numerical experiments. 


1 Introduction 

Partial differential equations with data varying on multiple scales in space and time, so called 
multiscale problems, appear in many areas of science and engineering. Two of the most prominent 
examples are composite materials and flow in a porous medium. Standard numerical techniques 
may perform arbitrarily bad for multiscale problems, since the convergence rely on smoothness 
of the solution [3]. Also adaptive techniques [26], where local singularities are resolved by lo¬ 
cal mesh refinement, fail for multiscale problems since the roughness of the data is often not 
localized in space. As a remedy against this issue generalized finite element methods and other 
related multiscale techniques have been developed [niiijiiniiTiisiiiniiiniiiniiii]- So far these 
techniques have focused on multiscale coefficients in general and multiscale diffusion in partic¬ 
ular. Significantly less work has been directed towards handling a computational domain with 
multiscale boundary. However, in many applications including voids and cracks in materials and 
rough surfaces, multiscale behavior emanates from the complex geometry of the computational 
domain. Furthermore, the classical multiscale methods mentioned above aim at, in different 
ways, upscaling the multiscale data to a coarse scale where it is possible to solve the equation to 
a reasonable computational cost. However, these techniques typically assume that the represen¬ 
tation of the computational domain is the same on the coarse and fine scale. In practice this is 
very difficult to achieve unless the computational domain has a simple shape, which is not the 
case in many practical applications. 

In this paper we design a multiscale method for problems with complex computational do¬ 
main. In order to simplify the presentation we neglect multiscale coefficients in the analysis even 
though the methodology directly extends to this situation. The proposed algorithm is based on 
the localized orthogonal decomposition (LOD) technique presented in [5D] and further developed 
in [71I1III1111]. In LOD both test and trail spaces are decomposed into a multiscale space and 
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a reminder space that are orthogonal with respect to the scalar product induced by the bilinear 
form of the problem considered. In this paper we propose and analyze how the modified multi¬ 
scale basis functions can be blended with standard hnite element basis functions, allowing them 
to be used only close to the complex boundary. We prove optimal convergence and show that 
the condition number of the resulting coarse system of equations scales at an optimal rate with 
the mesh sizes. The gain of this approach is that the global solution is computed on the coarse 
scale, with the accuracy of the fine scale. Also, all localized fine scale computations needed to 
enrich the standard finite element basis are localized and can thus be done in parallel. 

Other work on multiscale methods for complex/rough domains are |17LI18j . which are based 
on the multiscale hnite element method (MsFEM). However, the analysis in [TTIHSI is limited to 
periodic data. An explicit methods to handle complex geometry is the cut hnite element method 
(CutFEM) [5] which use a robust Nitsche’s formulation to weakly enforce the boundary/interface 
conditions. See also |12j for an other hctitious domain methods. An other explicit approach 
is the composite hnite element method [13] which constructs a coarse basis that is htted to 
the boundary. The hctitious domain and composite hnite element approaches are useful in 
e.g. multigrid methods since they construct a coarse representation on domains with hne scale 
features, see also m- The technique proposed in this paper is more related to the extended 
hnite element method (XFEM) [11] where the polynomial approximation space is enriched with 
non-polynomial functions. The method can be used as a discrete alternative to XFEM, that can 
be useful e.g. when the nature of the singularities are unknown. 

The outline of the paper is as follows. In Section we present the model problem and 
introduce some notation. In Section]^ we formulate a multiscale method for problems where the 
mesh does not resolve the boundary. In Section [^ we analyse the proposed method in several 
different steps and hnally prove a bound of the error in energy norm, which shows that the error is 
of the same order as the standard finite element method on the coarse mesh for smooth problem. 
In Section we shortly describe the implementation of the method and prove a bound of the 
condition number of the stiffness matrix. In Section]^ we present some numerical experiment to 
verify the convergence rate and conditioning of the proposed method. Finally in the Appendix 
we prove a technical Lemma. 


2 Preliminaries 

In this section we present a model problem, introduce some notation, and define a reference finite 
element discretization of the model problem. 

2.1 Model problem 

We consider the Poisson equation in a bounded polygonal/polyhedral domain C for d = 2, 3, 
with a complex/hne scale boundary dVl = Fd U F^j. That is, we consider 

—Au = f in n, 

iy-\7u + KU = 0 on Ffl;, (2-1) 

u = 0 on Fd , 

where ly is the exterior unit normal of 9H, 0 < k S M, and / S For simplicity we assume 

that, if /t = 0 then F^i yf 0 to guarantee existence and uniqueness of the solution u. The weak 
form of the partial differential equation reads: find u & V := {v H^{Vl) \ 'cItd = 0} such that 

a{u,v) := / VM-Vuda:-|- / KuvdS= / fvdx=:F{v), 

Jn Jtr Jq 


( 2 . 2 ) 
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for all V €V. Throughout the paper we use standard notation for Sobolev spaces [T]. We denote 
the local energy and L^-norm in a subset w C by 

f \ 1/2 

I'clIL = I / |VMpda;+ / Ku^dS] , (2.3) 


and 


respectively. Moreover, if w = we omit the subscript, |||u||| := |||u|||a and \H ■■= IMIn- 


'FrOBuj 
1/2 

'dx' 


(2.4) 


2.2 The reference finite element method 

We embed the domain in a polygonal domain fig equipped with a quasi-uniform and shape 
regular mesh Tf/,o, i-e., fl C tlo and tlo = J^tgTh o'^' mesh of Th,o 

consisting of elements that are cut or covered by the physical domain tl, i.e. 


TH = {TGTH^o\Tnn^(li}. 

The finite element space on Th is defined by 

Vh = {vG C°{n) I VT e Th, v\t G ViiT)}, 


(2.5) 


( 2 . 6 ) 


where 7^i(T) is the space of polynomials of total degree < 1 on T. We have Vh = span{(p 2 ,} 2 :gAt, 
where Af is the set of all nodes in the mesh Th and (p^ is the linear nodal basis function associated 
with node x G N. 

The space Vh will not be sufficiently fine to represent the boundary data. We therefore enrich 
the space Vh close to the boundary dH,. In order to construct the enrichment we define L-layer 
patches around the boundary recursively as follows 


w? :=mt{{T gTh \ Tr\n^T)r\n), 
Jr ■= int ((T G Tff I T n Jf^ ^ 0) n fl) 


for £ = 1,..., L. 


(2.7) 


Note that Wp is the set of all elements which are cut by the domain bou ndar y dVl. An illustration 
of n, Tnfi, Wp, and Wp are given in Figurej^ We will later see, in Lemma 4.8 that the appropriate 
number of layers is determined by the decay of the iL^-norm of the exact solution away from the 
boundary. 

Let 7)i be a fine mesh defined on J, obtained by refining the coarse mesh in J. On the 
interior part of the boundary dj \ dV, we allow hanging nodes. We define Vh{J) = {u G 
C°(a;p) I VT G Th, v\t G Vi{T)} and a reference finite element space by 

V^=(Vj/+V„(4))ni7^,(L!), (2.8) 

which consists of the standard finite element space enriched with a locally finer finite element 
space in J. We assume that the space is fine enough to resolve the boundary, i.e., we assume 
that the boundary dVl is exactly represented by the fine mesh Th- 

The finite element method posed in the enriched space reads: find Uh G such that 

a{uh, v) = F{v) for all v G V[. (2.9) 


We call the solution to (2.9) the reference solution. We have the following a priori error estimate 
|||w — < C(-ff+ h 


( 2 . 10 ) 


where 1 < s < 2 depends on the regularity of u in Wp For a proof of (|2.10[) see Section 4.3 
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Figure 1: The computational domain fl embedded in the mesh Th,o- The dark grey area is Wp, 
and the union of the dark and light grey area is cup. 


3 The multiscale method 

In the multiscale method we want to construct a coarse scale approximation of Uh, which can be 
computed at a low cost. We present the method in two steps: 

• First, we construct a global multiscale method using a corrected coarse basis which takes 
the fine scale variation of the boundary into account. 

• Then, we construct a localized multiscale method where the corrected basis is computed 
on localized patches. 

3.1 Global multiscale method 

For each x G JV (the set of free nodes) we define a L-layer nodal patch recursively by letting 

ui =: int {{f gTh , for £= 1,...,L. 

We consider a projective Clement interpolation operator defined by 

^HV = ^ {Pxv){x)(px, (3.2) 

xejGi 

where A/} is the index set of all interior nodes in 12 and is a local L^-projection defined by: 
find PxV G {v G Vh \ supp(u) C a;° 0} such that 

{PxV,w)^o = {v,w)^o for all w G {v G Vh \ supp(u) na;° 7 ^ 0}. (3.3) 

This is not the same operator as proposed in the original LOD paper [5D]. We choose the 
projective Clement interpolation operator since projective property simplifies the analysis and 
that it is slightly more stable for multiscale problems [ 1 ]. 
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Using the interpolation operator we split the space into the range and the kernel of the 
interpolation operator, i.e., Vh = and V-^ = (1 — Xh)V\. Since the space Vh does not 

have the sufficient approximation properties we use the same idea as in LOD and construct an 
orthogonal splitting with respect to the bilinear form. We define the corrected coarse space as 


= (1 + Q)XhVI 


(3.4) 


where the operator Q is defined as follows: given vh G Vh find Q(vh) G {n G | 'cItd 
such that 

a{Q{vH),w) = —a{vH,w) for all v G {v G \ tc|rD = 0}. 


-Vh} 

(3.5) 


Note that Vh V^ but V^ Gl because the correctors are solved with boundary conditions 
that compensates for the nonconformity of the space Vh and also that (3(i’i/)|a\(^*; = 0 since V^ 
only has support in Wp. From (3.5) we have the orthogonality a{Vjj,VV) — 0 and we can write 


the reference space as the direct sum = V]p 0 a V-^, where the orthogonality is with respect 
to the bilinear form a. 


The multiscale method posed in the space reads: find uh G VD such that 


a{uH,v) = F{v) for all v G V^. 


(3.6) 


Remark 3.1. Note that, even in the g loba l multiscale method the support of the correctors are 
only in Wp which we defined in Section 


2.2 


3.2 Localized multiscale method 

Finally, we further localize the computation of the corrected basis functions to nodal patches on 
Wp. Using linearity of the operator Q{vh), we obtain 

Qivn) = ^ vh{x)Q{ipx). (3.7) 

We denote the localized corrector by 

Q^{vh)= vh{x)Q’^{‘p^). (3.8) 

xeAfi 

where Q^^ifix) is the localization of Q{(px) computed on an L-layer patch. The local correctors 
are computed as follows: given x G Ni find Qxi^x) G {n G | = 0 and n|rD = —^x} 

such that, 

a{Q}^{(fx),w) =-a{ipx,w) for all G {n G | = 0 and ?;|ro = 0}. (3.9) 

The localized multiscale method reads: find G V^’^ := span{(/? 3 ; + Qx(^x)}xeAfi such that 

a{u^,v) = F{v) for all n G (3.10) 

r Z/ 

The space Vpp’ has the same dimension as the coarse space Vh but the basis functions have 
slightly larger support. The multiscale solution ujj G Vpp’ has better approximation properties 
than the standard Hnite element solution on the same mesh satisfy 

\\\uh-u}j\\\<CiH, (3.11) 

where is the mesh size, Ci is a constant independent of the fine scale features of the boundary 
F, and L = [(72 log 2 (iF“^)] for a constant C 2 (= 1-5 in the numercal experminents). A proof is 
given in Theorem |4.12| in the coming section. 
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4 Error estimates 


In this section we derive our main error estimates. First we present some technical tools needed 
to prove the main result: 

• We present an explicit way to compute an upper bound for Poincare-Friedrichs constants 
on complex domains. 

• We prove approximation properties of the interpolation operator on these domains. 

The main result is obtained in four steps: 

• We bound the difference between the analytic solution and the reference finite element 
solution |||it — Wft|||. 

• We bound the difference between the reference finite element solution and the ideal multi¬ 
scale method 

• We bound the difference between a function v S Vh modified by the global corrector and 
localized corrector |||Q(u) — ( 5 ^(?;)|||. 

• Together these properties are used to estimate the error between the analytic solution and 
the localized multiscale approximation. 

Furthermore, let o 6 abbreviate the inequality a < Cb where C is any generic positive constant 
independent on the domain and of the the coarse and fine mesh sizes Ff, h. 


4.1 Poincare-Friedrichs inequality on complex domains 

A crucial part of the proof is a Poincare-Friedrichs inequality with a constant of moderate size. 
The inequality reads: for all u S it holds, 


inf |ju — c||cj < C(a;)diam(w)||Vu||,, 


where the optimal constant is 


c = 


u dx. 


Following [53], we consider inequalities of the following type: for all u G 

||m - A.y(M)||,j < C'(w)diam(w)|| Vm||;^, 
where 7 C duj is a (d — l)-dimensional manifold and 

X^{u) = -—- [ udS. 

I7I Jj 


(4.1) 

(4.2) 

(4.3) 

(4.4) 


We introduce the notation Cpp = and refer to Cpp as the Poincare-Friedrichs constant, 

which depends on the domain ui but not on its diameter. 

A direct consequence of (4.3) is the following inequality 


IklL < C'pFdiam(w)||VM||^ -f diam(w)^/^||M||..y, 


( 4 . 5 ) 
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which holds if diam(a;)‘^ ^ ;< I 7 I, i.e., the average is taken over a large enough manifold 7 C du). 
A short proof is given by 


||m|U < IIm - A^(m)||^ + ||A^(m)||^ 

< C'pFdiam(w)|jVu||(^ + |A^(m)| • 

< C'pFdiam(w)|jVu||(^ + ||A^(w)||^|7|“^/^|w|^/^ 

< CpFdiam(a;)||Vu||tj + diam(a;)^/^||A-y(w)||-y. 


Furthermore, from [5^ we have the bound CpF < 1 for the Poincare constant on a d-dimensional 
simplex where 7 is one of the facets. 

Next we will review some results given in | 2 H] applied to domains with complex boundary. 
In |23j the notion of quasi-monotone paths is use to prove weighted Poincare-Friedrichs type 
inequalities using average on (d — l)-dimensional manifolds 7 C w. These results have also been 
discussed for perforated domains in [ 1 ]. 


Definition 4.1. For simplicity we assume that w is a polygonal domain that is subdivided into a 
quasi-uniform partition of simplices r = We call the region = {T■ -VJT 

a path, if 7^, and share a common (d — l)-dimensional manifold. We will call := s 

the length of the path and 77 = max 7 ’GT{diam(T)}. 

Lemma 4.2. Given t from Definition\4. 1\ and the index set = {£ : dTi H 7 0} it holds 


C/pF 


< 


^max^max^ 

\l\H^ 


d +1 


(4.7) 


where Smax = niax(sfcj) is the length of the longest path and rmax = nrax^gx \{{s,k) & I x J \ 
Ti G Pk,j}\ is the maximum number of times the paths intersect. 

Proof. See [53]. □ 

We will now use Lemma [42] to show some cases when Cpf can be bounded independent of 
the complex/fine scale boundary dfl. 


Fractal domain. We consider the fractal shaped domain given in Figure]^ First we compute 
Smax- The number of Ti on 7 is then proportional to 2^, where k is the total number of uniform 
refinements of the domain and we bound the maximum path length as 

Smax -- ^ < 2 • 2^ (4.8) 

i=0 ^ 

i.e., the maximum length of a path is proportional to 2^. Next we compute the maximum number 
of times a simplex is in a path, rmax- First we show how many times the elements besides 7 are 
in a path and then we show that this number is larger than on any other 7 ^, see Figure]^ The 
number of paths on each Ti is the total number of elements in the domain. On 7 we get 


r 


7 


fc k ,(0k]2 

'^na{i)eaii) = 



< 4•4^ 


(4.9) 


where «.□(*) is the number sub domains with index i and e\z\{i) is the number of elements inside 
a single sub domain with index i. Next we show that there is no T in other parts of the domain 






Figure 2: A fractal domain that has a bounded Poincare-Friedrichs constant. 


where the number of paths is proportional to something with a stronger dependence on n than 
r^. We obtain, 


r 




23 


k 

^nn(i)en(i) < 
i^j 


' 75 


(4.10) 


where 2-1 comes from that 13 2 ^ times smaller than 7 and 1/3-1 since the area of the domain 
affecting boundary r^j is less than 3-^ of the total area. This proves that rmax choosing 

structured paths in the interior of the squares. To finish the argument we note that H/rj = 2^ 
and I 7 I = H, and we obtain 


OpF 


< 


'^max^niax^ 

\l\H^ 


d+1 


< 1. 


(4.11) 


Saw tooth domain. An other example of a complex geometry is the saw domain given in 
Figure Let the width of the saw teeth be 7 = 2~^. A mesh constructed using 2^ uniform 



Figure 3: Saw tooth domain that has a bounded Poincare-Friedrichs constant. Here rj is the 
width of one of the saw teeth. 

refinements are needed to resolve the saw teeth. It is clear that Smax ~ 2^ and choosing the 
structural paths we have that r„iax ~ (2^)^. Again we have that Cpp ;< 1 as long the length of 
the saw teeth are fixed. 
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An example of a domain with a non-bounded Poincare-Friedrichs constant is e.g. a dumbbell 
domain. 

4.2 Estimation of the interpolation error 

In this section we compute the interpolation error for a class of fine scale functions needed in the 
analysis. For each T S Th we define an L-layer element patch recursively as 


uj^-=Tr\ n, 

Jt ■= int ((T e Tff I f n 7 ^ 0) n n) , for £ = l,..., L. 


(4.12) 


Lemma 4.3. The projective Clement type operator inherits the local approximation and stability 
properties for all interior elements, i.e., for all v G 


H-\v-Ihv)\\t + WVIhvWt < llVuLi , 


(4.13) 


holds for all interior elements T G Th- 


Proof. It follows directly from the standard proof [1] since Ti is a partition of unity on 

interior elements. □ 

The trace of a function v GV^ is “small” since the function v is in the kernel of an averaging 

operator, 

Vf = {vGVI\ Ihv = 0}. (4.14) 

We formulate this more precisely in the following Lemma. 

Lemma 4.4. Given an interior element T G Th let 7 C dT be one of its faces. Then 

\\vh < (4.15) 

holds for all v GVT 

Proof. For an interior element T the standard approximation property of the Clement type 
interpolation operator holds, i.e.. 


||t = \\v-Xhv\\t < H\\Tv\\^i. 


since v G . Using a trace inequality and (4.16) we obtain 


\\v\\i < H-^\\vrT+m^v\\T < 


where T is an interior element. 


(4.16) 

(4.17) 

□ 


We will now make an assumption which is a sufficient condition to prove the main results of 
the paper and which will also simplifies the analysis. 

Assumption 4.5. All elements in S' S Th share a vertex with an interior element, i.e., an 
element T G Th such that T HTl = T. 

Lemma 4.6. Let T be an element that is cut by the boundary dfl. Under Assumption \4.t^ the 
following Poincare-Friedrichs type inequality holds 


It for all v G . 


(4.18) 



10 




Figure 4: Admissible (left) and non-admissible (right) mesh according to Assumption 4.5 The 
irregular curve is where the elements are cut by the outer boundary. 


Proof. Let T be an element which share the vertex x with T. Using (4.5) we have that 


MIt < IkLo < H\\Wv\U + < H\\Wv\U < H\\WvL., 


(4.19) 


holds, since diam(a;°) < H, diam(w~) < H and C w- C wf, choosing 7 as an edge of T which 
contains the vertex x. □ 


Next we prove local approximation and stability properties for functions which are in some 
sense close to but outside V-^. 


Lemma 4.7. Letlu : L^(JIh) be the Clement interpolation operator defined by (3.2). If 

V = r]w, where rj and w satisfies 0 < 7 < 1, ||V77||z,oo < H~^, and w G then the following 
estimate holds 


\H-\v-lHv)\\T + \mHv\\T<\\yw\\ 


(4.20) 


for all T G Th- 


Proof. The local approximation and stability properties for an interior element follows directly 
from Lemma |4.6| together with 

I|V77^c||t < H-^\\w\\t + llVrullr < ||Vu;L^. (4.21) 

Next we investigate the local approximation and stability properties for elements on the bound¬ 
ary. Let T be an element cut by the boundary and T an interior element sharing vertex x with 
T. Then the L^-stability follows directly from the stability of the interior elements, i.e. 


||(P.u)(a:)||T = |(P.t^)(cc)|^||l||f < mv){x)\\\l\\f = ||(P.u)(x)||f 

<H^^^\\{P,v){x)\\^^^f^ < < WP^vWf 

< \\Pxv\\u>o < < IkL^- 


We obtain 


'-Ihv\\t < 


,1 < 


5 , iFllVr 


(4.23) 


since w GV^ using Lemma 4.4 L^-stability of the interpolation operator, and w GVC Similar 
argument yield 

\mHv\\T< H-^iP,v){x)\\T<\\^wL2, (4 24 ) 


where A/'t is all vertices in element T. 


□ 



















11 


4.3 Estimation of the error in the reference finite element solution 

We bound the error in the reference finite element solution Uh € V^. 


Lemma 4.8. Let u GV and Uh € be the solutions to (2.2) and (2.9) respectively. Then 


Iw—rt/illl < inf {\\H ^(u — vh) 
vneVn' 


+ lll«-^J^IIIn\^f-M+ inf 


\U-Vh\ 


(4.25) 


holds. 


Proof. We split LI into the different parts fl \ Wp, Wp \ Wp and Wp Since V[ C V, we have 
the best approximation result 


\U - Uh\ 


< 


for all w S VF. 


(4.26) 


Let ?7 G Vf/ be a cut off function, where = 0, ri\^k-i = 1, and || V? 7 ||pc=o(p) ^ H We 

construct w = vh + r-hr]{vh — vh) G where vh G Vh, Vh € Vh(wp) , and tt/j is the nodal 
interpolant onto the finite element space V/i and obtain 


|u — rcIlP = \\\u — vh'""^ 


Q\cl) 


k+\\\u-VH- TThVivh - + ll|n - 'I’hlWlk-l. 


(4.27) 


The first and third term are in the right form, see the statement of Lemma |4.8[ We turn to 
the second term. Using the fact that the nodal interpolant tt/j is iL^-stable for finite polynomial 
degrees (2 in our case) we obtain 


WWhViVh - VH)\\\\^,k-l < \\\vi.Vh - vh)\\\1 


I i-\ /c — 1 ^ Irl ) \\\ k\ fc — 1 

'a;p\u;r ~ in 'v ^' "cjp\a;p 

= \\^{vivn-vn))\\lk^^k-k 




(4.28) 


We have 


||V(r 7 (u;, 


< + Wv'^ivh - 

< H-^\\vh - + l|V(?;?, - VH)\\^ky^k-i 

< H-'^\\vh - u + U - VH\\^k\^k-l + \\V{Vh-U + U-VH)\\^ks^^k-l 


(4.29) 


rV 


<H + ||V(m-u//)||^ 




H 


- 1 | 






Furthermore 




< kH \\vh - + i^||V(^;;, - VH)\lk\^k-i. 

Taking the infimum and using that 

inf + ||V(u-?;?,)||^.\^.-i) 

vheVh \ r\r r\r/ 

since h < H and Vh C 14, concludes the proof. 


(4.30) 


(4.31) 


□ 


The analysis extends to a non-polygonal boundary if we assume that h is fine enough to 
approximate the boundary using interpolation onto piecewise affine functions, see e.g. |25] . 




4.4 Estimation of the error in the global multiscale method 

In this section we present and analyze the method with non-localized correctors. 
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Lemma 4.9. Let Uh G solve (2.9| and uh G solve (3.6), then 

\\\uh-UH\\\ < H\\f\\^k 

holds. 


(4.32) 


Proof. Any G V[ can be uniquely written as Uh = uh + u^ where uh G and G V^. This 
follows from the result from functional analysis, that if we have a projection V : Uh ^ ur onto 
a closed subspace, there is a unique split Uh — Vuh + (1 — 'P)uh. For 7^ = (1 + Q)Ih, we have 
VVj^ = Vr and 

= (1 + Q)Ir{1 + Q)Ir = (1 + Q)IrIr = (1 + Q)Ir = V. (4.33) 

We obtain 

|||Mf|||2 = a{u\u^) = = {f,u^-lRU^)L 2 (a,^) < H\\f\\^k\\\u^\\\, (4.34) 

which concludes the proof. □ 


4.5 Estimation of the error between global and localized correction 

The correctors fulfill the following decay property. 


Lemma 4.10 (Decay of correctors). For any x G Ni there exist a 0 < 7 < 1 such that the local 
corrector G V[(ipx) o,nd the global corrector Q{(px) G V^{i.px), which solves (3.9) and 

(3.5) respectively, fulfills the decay property 


IIKQ - Q^){vr)\\\^ < L^L(i-3)/3J ^ vl\\\Qdgx\\\\ (4.35) 

xGTV/ 


where [-J is the floor function which maps a real number to the largest smaller integer. 

Proof. See Appendix [A] □ 

The localized corrected basis functions fulfill the following stability property. 

Lemma 4.11. Under Assumption \4.t^ we have the stability 

\\\Fx + Q^iPxM < H-^\\(px\\, (4.36) 


for the corrected basis function given any x G JVj. 

Proof. First we will prove that there exist a (non-unique) function gx G such that 

idx — Px)\td =0 IIl 5 a;lII ^ Ill¥’a:llI for all s S Nr Given any node x define w\t = Px — Fx and 
■u^InXT = 0 where T is an interior element. The function w have to fulfill the following restriction 

Irw = iRipx = Fx, (4.37) 


which is equivalent to 


{Pyw){y) — Sxy, 


( 4 . 38 ) 
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where 6xy is the Kronecker delta function. In order to construct w we perform two a uniform 
refinements in 2D. A similar construction is possible in 3D using two red-green refinements but 
we restrict the discussion to 2D for simplicity. Then we have three free nodes in T for a function 
that is zero on the boundary dT. We write w as w = '^jVj where (pj are the Pi Lagrange 
basis function associated with the three interior nodes. We can determine w by letting it fulfill 

d+l 

'y ' ^ji^Vi^j)iyi) = ^x,yi- (4.39) 

i,i=l 


The value {Py(pj){x) can be computed as 

{PyPj){y) = (4.40) 


where Mu /a is a local mass matrix computed on = 1 for index x and 0 otherwise, and 

n : Vh Vh/a is a projection from the finer space onto the coarse. On a quasi-uniform mesh 
Pyipj(i/) is independent of H. Therefore, aj is independent of PI and there exist a constant such 
that 


Ikii < c\\ipx\\- 


(4.41) 


This yields 


\w\ 


< 4:H 


-11 


\w\ 




(4.42) 


Using the triangle inequality we have 


ll|ga;lll < lllv’xlll + infill < iL '^WPxW- (4.43) 

Next we consider the problem: find Qo G such that 

a{Qo,w) = a{(px - gx,w) w G V'^(w^), (4.44) 


where gx satisfies gx G {px-gx)\TD = 0: and |||gx||| < lll</^x|||- It is clear that Q^{px) = 

Qo -f g. For the stability we obtain 

\\\px + Q^iPxM"^ < a{(px + Q^ipx),‘Px + Q^[px)) = a{ipx + Q^{px),px + Qo+g) 

= a{(px + Q^iPx),Px+ g) < lllv^x + Q^{Px)\\\{\\\Px\\\ + Ill5a:|||) (4.45) 

< WIPx + {PxM ■ Ilk:i:|||, 

which concludes the proof. □ 


4.6 Estimation of the error for the localized multiscale method 

The a priori results for the localized multiscale method reads. 


the bound 


Theorem 4.12. Let u G V solve (2.2| and G solve (3.6). Then under Assumption 4-5 


I u — 111 < inf 

vGVh 


H '^{u — v) 


inv 


k-i -|- m — u 


I n\:». 


+ inf \\\u-v\\\^k + \\HfLk+L<^/^H~^j^\\f\\, 
veVhioj^) ^ 


(4.46) 


holds for k > 2. 
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• r Z/ 

Proof. Since Vjf C V we have the best approximation result 

r rZ/ 

Up ~^ Up ~'*^-^111 ior a\\ VH (z Vif . 

We obtain 

|||m - u^lll < |||m - Mhlll + - M/i-lll + |||w/i- - r)_H-|||, 


(4.47) 


(4.48) 


using the triangle inequality. The first and second term is bounded using Lemma 4.8 and |4.9| 
For the third term we choose vh = uh + Q{uh) which gives 


\UH - vh\\\^ = IIIQ(wff) - Q^{uh)\\\^ < ^ vl\\\Q(i), 

xGAfi 


2 

^ I I I 5 


(4.49) 


using Lemma |4.10| Lemma |4 .11 1 now gives 


\UH - VhIW^ < H ^ UHix)‘^\\Px\\‘^ < L‘^^‘^^\\uh\\ 

x^Mi 


(4.50) 


where also a Poincare-Friedrich inequality has been used. Combining (4.48) and (4.50) concludes 
the proof. □ 


5 Implementation and conditioning 

In this section we will shortly discuss how to implement the proposed method and analyze the 
conditioning of the matrices. For a more detailed discussion on implementation of LOD see m- 


5.1 Implementation 

To compute Q^(}Px) in ( |3.9[ ) we impose the extra condition = 0 using Lagrangian multipliers. 
Let Ux and Nx be the number of fine and coarse degrees of freedom in the patch a;°. Let Mx 
and Kx denote the local mass and stiffness matrix on ui^ satisfying 

{v,w )^0 vFMxV, (5.1) 


and 


o(u, w)|^o w'^KxV, 


(5.2) 


where v,w £ V/t|^o and w,v G are the nodal values of v,w. We also define the projection 
matrix Ha; : {z; G V//(w°) | supp(u) n a;° 7 ^ 0 } —)■ V/t(a;°) of size {ux x Nx) which project 
a coarse function onto the fine mesh and a Kronecker delta vector of size Nx x 1 as (5a; = 
(0,..., 0,1,0,..., 0) where the 1 is in node x. We obtain 


Pxv{x)=0 Xlv = Sx(nlMxUx) ^UIMxV = 0. 


(5.3) 


Set A = [Aya, Aya, • ■ ■, then (3.9) is equivalent to solving the linear system 


'Kx A 

Q^iPx) 


^x 

0 

. . 


0 


(5.4) 


where Q^{(px) G contains the nodal values of Q^{(px) and /r is a Lagrange multiplier. For 
each coarse node x in wf we need to invert n^Ma;na; to assemble (5.4), however since the size is 
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only (Nj: X Nx) it is cheap. The coarse scale stiffness matrix K in (3.10) is given element wise 

by 

= a{ipj + Q^{(pj), (fii + Q^{ipi)). (5.5) 

We save the nodal values of the corrected basis 


$ = 


+ Q + Q iV’y2 ) ■ • ■ i V^yjv + Q {VyN ) 


(5.6) 


Given the fine scale stiffness matrix K and the collection of corrected basis functions $, we can 
compute the coarse stiffness matrix as 

K = (5.7) 


and the linear system (3.10) is computed as 


= b, 


(5.8) 


where is the nodal values of and b correspond to the right hand side by^ = (/, ipy^ + 
Q^i^Pvi))- However, the fine stiffness matrix does not have to be assembled globally. If a Petrov- 
Galerkin formulation is used, further savings can be made [5]. 


5.2 Conditioning 

The Euclidean matrix norm is defined as 


N = sup 


\Av\ 


N 


(5.9) 


where {v,w) = ^^d \v\n = \/{v,v). 

Theorem 5.1. The bound 

k=\\K\\n\\K-^\\n <H-^, (5.10) 

on the condition number k holds. 

Proof. To prove the bound of the condition number we use the following three properties. 

1. An inverse type inequality for the modified basis functions. We have 

\\\T^ + QiTi)\\\^ (5.11) 

from Lemma 14.111 


2. A Poincare-Friedrich type inequality on the full domain, see Section 4.1 

3. An equivalence between the Euclidean norm and the L^-norm. We have that 

< E 


(5.12) 
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and 


N 


N 


N 


\vfN = ^ ^ = H-^W^Hvf < H-%f. (5.13) 

i^l i^l i^l 

holds, hence \v\n ~ H~'^/'^\\v\\. 

We have 

\Kv\„= sup ™ IIHMIHII ^ „u-s 


= sup 


= sup 


P|iv O^wGR^ 

using property 1) and 3). Also 

\K-\\jf = H-'^WK-^f < CpfH-^^IWK-^IW^ 

< H-‘^{KK-\,K-^v) < H-‘‘\v\n ■ \K-^v\n, 

using property 2) and 3). The proof is concluded by taking the supremum over v. 


(5.14) 


(5.15) 


□ 


6 Numerical experiments 

In the following section we present some numerical experiments that verifies our analytical results. 
In all the experiments we fix the right hand side to / = 1 and the localization parameter 
L= ri.51og2(i?-i)l. 

6.1 Accuracy on fractal shaped domain 

We consider the domain in Figure We use homogeneous Dirichlet boundary conditions on the 
most left, down, and right hand side boundaries and Robin boundary condition with k = 10 
on the rest, see Figure We compute localized correctors in the full domain. The reference 
solution is computed using h = 2“®. As seen in Figure a even higher convergence rate than 
linear convergence to the reference solution is observed. We also see that the correct scaling of 
the condition number with respect to the coarse mesh parameter H. 


6.2 Locally added correctors around singularities 

Let us consider two different domains, a domain with a re-entrant corner ([0,1] x [0,1])\([0.5,1] x 
[0,0.5]) and a slit-domain ([0,1] x [0,1])\([0.5,0.5] x [0,0.5]) with homogeneous Dirichlet boundary 
condition. We only compute correctors in a vicinity of the singularities, a nd d enote this domain 
The size of uj^ can be determined by the size of \u\h 2 (uj^) (see Section 2.2), i.e., how fast the 
singularity decay away from the boundary. We choose as in Figure 

As seen in Figure the correct convergence rates to the reference solutions and the correct 
scaling of the condition number are observed for both singularities. 


6.3 Locally add correctors around saw tooth boundary 

Let us consider a unit square where one of the boundaries are cut as a saw tooth and where 
correctors are only computed in a vicinity of the saw teeth shown in Figure ]^ On all the non 
saw tooth boundaries we use homogeneous Dirichlet boundary conditions. On the saw tooth 
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Figure 5: The fractal computational domain used in Section [6.1| We use homogeneous Dirichlet 
boundary conditions on the boundary segments marked by a thicker line and Robin boundary 
conditions elsewhere. 


boundary we test both homogeneous Dirichlet and Neumann boundary condition. We observe 
the correct convergence (although the Neumann case is a bit more sensitive and might need 
slightly larger domain where correctors are computed) and scaling of the condition number, see 
Figure 

6.4 Accuracy and conditioning on cut domains 

Let us consider a domain with a re-entrant corner D = ([0,1] x [0,1])\([0.5,1] x [0,0.5]). We want 
to investigate how the location of the boundary relative to the coarse background mesh affects the 
accuracy in the approximation and the conditioning of the matrix. We fix the coarse H = 2~^ and 
fine h = mesh sizes and consider three different setups of boundary condition, homogeneous 
Dirichlet on the whole boundary, homogeneous Dirichlet on the cut elements and homogeneous 
Neumann otherwise, and homogeneous Neumann on the cut elements and homogeneous Dirichlet 
otherwise. We will cut the coarse mesh in two different ways, 1) with a horizontal cut and 2) 
a circular cut around the reentered corner, see Figure [m The errors are measured in energy 
norm. The results are presented in Table and We conclude that neither the error nor the 
conditioning are sensitive to how the domain is cut. 
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Mesh size H 


Figure 6: The convergence rate to the reference solution in relative energy norm (left) and the 
scaling of the condition number (right) for the fractal domain in Figure]^ 




Figure 7: The marked area w* is where the finite element space is enriched, for domain with a 
re-entrant corner (left) and domain with a slit (right). 


r cut 

on \ 

erel(a) 

erel(6) 

Crel (c) 

D 

D 

0.059 

0.057 

0.056 

D 

N 

0.018 

0.019 

0.020 

N 

D 

0.063 

0.055 

0.053 



cond(a) 

cond(5) 

cond(c) 

D 

D 

9.85 

10.10 

13.63 

D 

N 

299.75 

282.26 

353.27 

N 

D 

10.537 

10.79 

11.47 


Table 1: For cut 1 we have L n [0,1 — r] x [0,1], where r = {/i, 0.5i7, 77 — h} in a), 6), and c), 
respectively. The errors measured in relative energy norm and condition number of the coarse 
stiffness matrix are presented. We try the different boundary conditions, D (Dirichlet) and N 
(Neumann), on the boundary segement Fcut, which cuts the elements. 


A Proofs 


In the appendix we present the proof of Lemma 4.10 
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Figure 8: The convergence rate to the reference solution in relative energy norm (left) and scaling 
of the condition number (right) for the re-entrant corner and the slit domain. 




Figure 9: We consider the saw tooth domain (left). The marked area (right) is where the finite 
element space is enriched. 


r cut 

dn \ 

erel(a) 

erel(6) 

^rel (c) 

D 

D 

0.060 

0.064 

0.073 

D 

N 

0.0205 

0.035 

0.048 

N 

D 

0.060 

0.057 

0.059 



cond(a) 

cond(6) 

cond(c) 

D 

D 

9.80 

9.23 

7.03 

D 

N 

246.99 

107.52 

59.67 

N 

D 

9.90 

11.44 

12.16 


Table 2: For cut 2 we have L\B{xo, r) for a ball B centered in xq = (0.5,0.5) and with radius r, 
where r = {/i, 0.5H, H} in a), b), and c), respectively. The errors measured in energy norm and 
condition number of the coarse stiffness matrix are presented. We try the different boundary 
conditions, D (Dirichlet) and N (Neumann), on the boundary segement Tcut^ which cuts the 
elements. 


Proof of Lemma \4-l(A We will make frequent use of the cut off function rj^ which satisfy 
jyfc-Ffc = 0 in = 1 in \ and ||< H-^. Given vh G Vh with 

nodal values Vx we define define e = {Q — Q^){vh) and e = where = {Q — Qx)vx<Px, 
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Figure 10: The convergence rate to the reference solution in relative energy norm (left) and 
the scaling of the condition number (right) for the saw tooth shaped boundary using different 
boundary condition on the saw tooths. 


we obtain 


|e||P < a(e, e) = a a{ex,e-Vx) 

xeN'i x^Mi 


(A.l) 


for Vx which satisfy 


a{ex,Vx)=0. (A.2) 

We choose = ry^+^^+ig _ 2 ^^l+ 2 ,l+i^ _|_ .sphere Xugx = 0, gx\dn = and 

|||5a:||| < ll |3ig7? )(^^A+ie|||. It is possible to construct gx with the given properties, see the proof 
of Lemma 4.11 Equation ( A.2[ ) now follows since Vx G and supp({j 2 ,) n supp((5,(r(u)) = 0. For 
each X € Mi, we obtain 

a{ex,e - Vx) < |||ea:||| • |||e - UxlH (A.3) 


where we split 


|e — Uxlll < ll|e — 




e-Vx 


The first term in dAAt can be bounded as 

111(1 - < |jv(i _ + 11^(1 _ 


')e| 


awl'+^nr 


(A.4) 


(A.5) 


< e 


using the product rule, inverse estimates, and interpolation estimates. The second term in (A.4) 
can be bounded as 


7 ^+^A+ig _ 73^|||2 ^ ^ ||^i^(^ie)||^ < |||e|||^+4 (A.6) 


where we used 


lVIj^(r?^+2A+ig)||2 ^ ||v(I^r?^+2.^+ie)||2,+3\^. < 


(A.7) 
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Figure 11: A given background mesh that is cut in two different ways 1) and 2) with different 
size of the cut a), b), and c). 


and 


X I * 33 '-‘^X 

(A.8) 

< kH ^||e — Iije|Pi+3 < KF^|jVe|Pi+4 < |||e|||^i,+4 


which follows from Lemma 4.7 and that k < H Hence, combining (A.l), (A.4), (A.5), and 

/ \ 1/2 


(A.6) we obtain 


lelll" < ^ me, 

xGAd 




<id E III'- 

\xeAfi 


(A.9) 


(A.IO) 


that is 

xeMi 

We use ic, = (1 - ri^~^’^~‘^)Q{vx(px) + ^H'nx~^’^~^Qivx'^x) + gx & V^uj^) where we use the 
same construction of Qx as for gx and obtain 




= a{ex, r]x ’ Q{vx^x) - ’ Qivx^x) + 9x) 


< 


L-l,L-2 


"XlWWW'lx 


Q{vxipx)\\\ + \\\IhVx "^Qivxg}. 


and hence 


^ ll|6a;||| ■ 11 IQ(l'a;V2a;) 11 

lllexlll E 111 |Q(‘/2a:) 11 


(A.ll) 


(A.12) 
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We ODtam 


dS” 


(A.13) 


Next we construct a recursive scheme which will be used to show the decay. 

< f r]^''^~^VQ(pa;'VQipxdx+ f dS” 

j Q, r fi 

= [ 'VQ(p^'V{ri!^’^~^Qip^)dx+ I K,Qifx{jl^^~^Qv>x) 

Jn xtr 

— / Q‘Px^Q‘Px^'i1t^~^ dx 

Jn 

= a{Qtpx,'n^’"~^Qifx) - I Q(Px'^Q‘Px'^r]T''~^ 

Jn 

in (A.13) can be bounded as 

'^QVx) = a{Qipx,Vx''"~^QVx -'^hVx''"~^QVx - 9x) + a{Q(px,lHrix’'^~^Q‘Px + Qx) 
= a{Q(px;,lHVx’’'~^Q‘Px + Jx) < \ \\Q(Px\\\^k\^k-l\\\lHVx’'~'^Q(Px\\\ 

^ \\\QVx\\\^k\^k-l\\\Qipx\\\\^k + l^^k-2\\\QiPx\\\ljk + l\^^k-2 

where again we construct such that |||ga;||| < \\\IhVx’^~^Q‘Px\\\ and use Lemma 
second term is bounded as 

/ Qipx"^Q^x"^Vx’''~^ dx < - InQxvW^kx^i^-AWQy^xlW^kx^i^-i 

Jn " (A.15) 


The first term 
a{Qv’x,Vx^~ 


4.7 The 




~ iii'vr^iMo;^ \u;^ m 

Combining (A.13), (A.14), and ( A.15| ) we obtain 


t-2 < |||(5l^x|||^fc + l^^fc-2 


11 11 ln\a;J — 11I 


^ :) 


and hence 




I I IV r-^ / I Mi2\cL)~ — /Ml^ 

where 7 = • Using ( |A.17 ) recursively we obtain 


,fc -2 , 


11 11 ln\w 

77 11 IQ^Pa: 11 la\,.,i.-3 




Combing (A.9), (A. 12 1 , and, (A. 18) we get 


which concludes the lemma. 


x^Mj 


(A.16) 

(A.17) 

(A.18) 

(A.19) 

□ 
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